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ON  THE  INITIATION  OF  DETONATION 
IN  A  ONE -DIMENSIONAL  SYSTEM 

This  Technical  Report  is  meant  only  for  internal  review, 
as  it  would  be  difficult  to  imagine  the  contents  are  not 
known  to  those  active  in  the  field.  The  purpose  in  writing 
it  is  to  provide  Institute  members  with  a  concise  discussion 
on  the  computation  of  the  flow  behind  a  shock  leading  to 
detonation  in  a  one-dimensional  system.  The  computation 
would  be  performed  using  the  method  of  characteristics,  and 
this  reduction  would  also  be  utilized  later  for  our  investiga¬ 
tion  of  the  stochastic  initiation  of  detonations.  Stochastic 
effects  might  be  introduced  either  as  white  noise  or  as 
variations  in  material  properties  such  as  so-called  "hot  spots." 
Mathematically  these  effects  would  appear  as  part  of  the 
"driving- terms"  in  the  characteristic  formulation.  Consequently 
the  purpose  of  the  present  Technical  Report  is  to  be  viewed 
primarily  as  a  supportive  document  for  our  research  effort  on 
the  stochastic  initiation  of  detonation. 

The  present  deterministic  system  is  supposed  to  consist 
of  material  whose  equation  of  state  resembles  that  of  a 
polytropic  gas.  We  assume  that  the  system  is  originally  shocked 
by  an  infinite  piston  moving  at  a  constant  velocity,  and  that 
this  shock  is  sufficient  to  start  certain  chemical  reactions. 

These  reactions  initiate  a  detonation  front  (shock  wave) ,  which 
in  turn  will  propagate  through  the  system  and  sustain  the  chemical 
reactions  in  a  zone  immediately  behind  the  front.  The  reaction 


is  assumed  to  transform  undetonated  substance  into  a  material 


with  similar  thermodynamic  properties,  plus  a  certain  amount 
of  heat  energy.  The  rate  at  which  heat  energy  is  released 
depends  on  the  local  thermodynamic,  and  progress  variables. 

At  this  point  we  introduce  the  notations  that  are  used 
in  this  work. 

x  -  distance  measured  from  an  original  point 
t  -  time  measured  from  the  initial  point 
u  -  particle  velocity 
p  -  material  density 
t  -  specific  volume,  t  *  ^ 

S  -  specific  entropy 
e  -  specific  internal  energy 
q  -  heat  energy 

i  -  specific  enthalpy,  i  -  e  +  pt 
c  -  velocity  of  sound 
p  -  pressure 
T  -  temperature 
B  -  adiabatic  constant 

g  -  specific  heat  at  constant  specific  volume 

T 

U  -  shock  velocity 

g  -  de/dt 

kl’  k2’  *1*  a2 ’  bl*  b2'  po’  UA’  Y-  are  constants 


H 


Given  Data 
UA  -  piston  velocity 
pQ  -  density  of  unshocked  gas 

Q  -  energy  of  formation  of  detonation  products 
kf,  -  Constants  in  the  rate  equation 
a^,  a2>  b^,  b2  -  constants  in  the  iteration  routines 
Y  -  equation  of  state  constant 

coo’ 

N,  -  number  of  points  on  each  C+  characteristic 


\/tu- 


We  also  introduce  the  following  notation  for  derivatives 


A  •  ”  -ir^  A  •  *  A- 
x'  3x*  Af  Tt* 


dA  -  SA 


3A 


Ht  =  Jt  +  u-5x  • 


1 .  The  Mathematical  Model: 

In  the  regions  of  space-time  in  which  shocks  are  absent  the 
flow  is  to  be  governed  by  the  following  conservation  equations. 


(mass)  Pt  +  upx  +  pux  -  0,  (1.1) 

(momentum)  u_  +  uu„  +  \  p„  *  0,  (1.2) 

(energy)  TS  -  q  ■  e  +  px  .  (1.3) 

Outside  of  the  reaction  zone  (3)  may  be  replaced  by  S  *  0; 


that  is  the  entropy  does  not  change  along  a  particle-path  or 
stream-line.  In  the  reaction  zone  it  is  convenient  to  express 
the  energy  balance  (1.3)  by  the  following. 
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q  Q6  =  e  +  pt,  _ (1.3’) 

and 

6  -  g(p  ,  p ,  S)  ,  _ (1.4) 


where  g(p,p,8)  is  a  function  of  the  density,  the  pressure,  and 
the  progress  variable  6  which  describes  the  chemical  kinetics 
of  the  material. 


The  questions  1.1,  1.2,  1.3',  1.4  are  a  system  of  first 
order  partial  differential  of  the  form, 


Lim  aij  +  bij  “Tt  +  ci  “ 


(1.5) 


(i  =  1.2, 3, 4)  ,  (j  =«  1,2, 3, 4)  , 

where  the  coefficients  a,  b,  c,  depend  on  x,  t,  $k.  A  necessary 
and  sufficient  condition  is  given  for  the  existence  of 
characteristic  directions  for  this  system  of  equations,  that 
is  conditions  (for  the  existence  of  curves  <7k,  such  that  a 
directional  derivative  may  be  formed  in  the  same  direction 

as  [1],  [2]  The  necessary  and  sufficient  conditions  for 
characteristic  direction  is  that  the  following  homogeneous 
equations  for  the  A±  be  compatible  [1],  [2] 


Xi(aijca  -  bijxo> 
Xi(aij^J  +  C;Lxa) 

Xi<bij*aj  +  V< P 


0, 

0, 

0, 


(1.6) 


5 


where  j  =  1,2, 3, 4.  If  one  sets  the  determinant  of  the  first 


four  equations 


aijc«  -  bijxa '  *  “■ 


we  may  obtain  an  algebraic  relation  for 
the  following  association 


12  3  4 

<t>  :  =  0,  f  :=  U,  <$  :=p,  f  :=  $, 


By  making 


and  formally  computing  this  determinant  one  obtains 

_  s2 


'•ij'o  -  bljX3l 


(ut  -x  )  0  „ 

--PTY---D--  ([Uto  ’V  ’  YtQ2  -  | 


'a  p 


(1.7) 


The  characteristics  are  then  given  by  the  equations 
l+,-:  uta  -  xa  “  *  ±cac-  _ 


(1.8) 


ZQ  :  (uta  “xa}  "  0 


(1.9) 


The  first  two  I+,  I_,  represent  the  paths  of  sound  waves ,  whereas 
the  second  two  (a  double  root),  Iq ,  are  degenerate  characteristics 
and  coincide  with  the  streamlines . 

The  remaining  equations  in  the  system  (1.6)  give  rise  to 
compatability  conditions  on  the  thermodynamic  variables  in 
order  that  (1.8)  and  (1.9)  be  characteristics.  Corresponding 
to  the  characteristics  I+  _  we  have  the  compatability  conditions, 

p 

n+,-  :  ±P%  +  -|  -  QgP  _ (1.10) 


(1.10) 
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and  corresponding  to  Iq ,  we  have 


II, 


Na  =  »V 


(1.11) 


Because  the  flow  is  non-adiabatic  in  the  reaction  zone 
(S  f  0)  one  can  not  introduce  Riemann  invariants  by  integrating 


2 .  On  the  Numerical  Computation  of  the  Characteristics: 

In  this  section  we  will  develop  a  method  by  which  the 
characteristics  can  be  approximated.  We  shall  assume  that 
two  neighboring  points  P2,  P3  are  given  on  an  I+  characteristic, 
and  that  P^  lies  on  the  I_  characteristic  through  ?2-  Further¬ 
more,  the  image  points*  in  the  (u,p) -plane,  ^ wiH 

also  be  assumed  known.  If  the  points  Pk  are  sufficiently 
close  the  characteristics,  which  join  them  may  in  general  be 
approximated  by  parabolic  arcs.  In  this  case,  we  may  write  that 
on  I+,  one  has 


1 

1 


f 

dx 

1  1 
1 

1 - 

u4+c4+ui+ci j 


(2.1) 


or  by  introducing  a  bar  and  double  subscript  for  the  average 
value  at  two  positions, 


*  For  each  point  P  s  (x,t) -plane  there  corresponds  an  image 
point 

J 

it  e  (u,p)-plane,  given  by  the  correspondence  u  :=u(x,t), 
p  :-p(x, t)  . 


FT 


- wtr  — *■ 


r 


i 
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I 


+ 


X4~XL 

t4-ti 


+  c 
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The  compatability  condition,  II+,  become 


u4~ul 

P4-P1 


+  2Q 
14 


fal 

icj 


t4‘t 

IP4-P1 


(2.3) 


On  I_,  and  II_  one  has  respectively 


and 


x4~x3  ;  -  -  , 

t^-t^  U43  "  c43 


u4~u3 

P4-P3 


U) 

IpcJ34 


t4"t3 
34IP4-P3 


(2.4) 


(2.5) 


In  order  to  obtain  a  first  estimate  of  u^,  p4,  S^,  we 
replace  the  secants  ,  by  the  tangent  lines  at  and 

that  is  by  the  lines  drawn  through  and  with  the 
respective  slopes 


and 


(2.6) 


(2.7) 


The  intersection  of  the  tangent  lines  yields  an  estimate 

JL 

for  the  point  i=(u4,  p4)  ,  which  we  use  to  estimate 

P4,  and  c4. 
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If  our  material  is  a  polvtropic  gas,  then  it  obeys  an 
equation  of  state  having  the  form, 


P 


(7-1) 07  exp 


rs-so' 


(2.8) 


where  Sq  an  appropriate  constant.  Since,  we  have  an  estimate 
for  P4,  knowledge  of  either  S  or  p  at  implies  knowledge  of 
the  other.  We  may  obtain  an  estimate  of  by  making  use  of 
the  degenerate  characteristics,  Iq  ,  IIq.  Along  IIq  we  have 


TS  = 

ot 


(7-1) P  C 


Sa  QSa 


(2.9) 


or 


,  I  Sn-S  I 

s»  *  eitpjmrj  st“  ’ 


(5.10) 


We  estimate  S4  *  S4*  by  assuming  the  streamline  passes 
through  and  P^,  and  then 


V~si  ;  hgi 
V'l  ofT  P 


■s  -S,7 

o  1 

T - 


2  v  -1 


(2.11) 


and  which  in  turn  allows  us  to  approximate  =  c^'v .  An  estimate 
for  c4  may  be  obtained  from 


f  3p 
I ’So 


*  ,  Q 

*  /  ,  v  ,  Y-l  I  4  0  ;  . 

■  7P4  “y(y-1)(p4  )t  exp  j — c - 1 

4  L  v  j 


(2.12) 


x4”xl 


point  P,  may  be  found  by  computing  the  ratios  - — — 
4  c4'cl 


x4-x3  *  ,u 
-  — g—  from  the  approximate  values  u4  ,  c4" . 
a4  2 
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We  recall,  that  in  order  to  approximate  at  we 
assumed  that  the  streamline  passed  through  P^,  P^.  To  improve 
this  estimate  we  may  compute  the  slope  of  Iq  at  P^,  and  obtain 
the  intersection  of  this  line  with  I+  or  I_ .  (See  Figure  2) . 
We  call  the  intersection  P^ ,  which  we  estimate  by 


5 

5 


* 

u4  =  u4 


(2.13) 


The  value  of  u,  S  at  P^  may  be  obtained  by  interpolation  using 
the  values  of  u,  S  at  P^,  ?2<  P3*  The  secant  of  the  IQ  through 
P^,  P^  then  has  the  slope 


(2.14) 


We  may  obtain  an  improved  approximation  of  S^,  now  using 
the  difference  expression. 


S4"S5  *  ^t4't5^ 


Cv(8) 


45 


exp 


Vs' 


v 


1 


(2.15) 


45 


where  averages  are  taken  from  values  obtained  at  P^  and  the 
first  estimates  at  P^.  With  this  value  of  we  can  obtain 
another  estimate  for  from  the  equation  of  state.  Using 
these  new  estimates  of  p^,  P^,  etc.  We  may  now  improve  our 
location  of  in  the  (u,p) -plane  by  means  of  the  difference 
formulae  for  secants  instead  of  the  tangent  approximation. 


10 


If  the  differences  between  first  and  second  estimates 

appears  too  large,  it  is  suggested  that  the  procedure  outlined 

above  be  repeated.  We  include  at  this  point  details  from  the 

scheme  to  compute  the  flow  behind  the  shock,  which  we  assumed 

is  strong.  Furthermore,  we  restrict  our  investigation  to  the 

f  -ko  ^ 

consideration  of  reaction  rates  of  the  form  g:=  exp  ^  - j  . 

Our  scheme  consists  of  three  routines  which  we  lable  (A) ,  (B) , 
(C) .  These  are  indicated  in  what  follows .  We  use 


(A)  for  general  points , 

(B)  for  piston  path  points, 

(C)  for  shock  points. 


We  intiate  our  calculations  at  point  A  (See  Figure  1) ,  and 
calculate 


U, 


^  UA-  PA  =  po  UA  UA’  PA 


1/2 


®A  "  0*  and  gA  *  kl  C'k2?A/PA  Then 
3  to  indicate  a  generic  variable,  i.e. 


let  Yoo  "  YA,  xoo 


Uoo  ^00 


Use  routine  (B)  to  get  the  solution  at 
(here  A-l,  oo-2,  10-3) 

Use  routine  (C)  to  get  the  solution  at 
N  parts  and  interpolate  for  the  values 


let  us  use  the  symbol 
Y  =  U,  u,p,o,c,8,g,  and 

point  10. 

point  IN.  Divide  into 
of  the  variables,  i.e., 


GENERAL  POINT 


ADVANCE 

ra 


TEST 


U-l)  U)  (i-l) 

al  a3  -a3  1  a2  a3 


a=x,  t ,  p,  p ,  3,  g 


ADVANCE  i 


NEXT  POINT 
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3 .  The  Characteristic  Equations  in  the  Steady  State. 

The  description  of  the  one-dimensional,  steady-state, 
detonation  of  an  ideal  gas  may  be  found  in  [4] .  In  this  section 
we  shall  use  their  solution  in  order  to  obtain  the  steady 
state  characteristics,  I  .  These  characteristics  will  be 
useful  to  us  in  solving  for  the  transcient  characteristics 
In  the  undetonated  region  the  specific  enthalpy  may  be 
represented  in  the  form  [3] 


i0  =  y±T  PoU0  +  Q  =  Vo  +  Q?  _ (3'L) 

right  after  the  shock  front  (but  before  any  reaction  occurs)  we 
represent  it  as 

i  =  yTT  PXU1  +  Q  =  Vl  +  Q*  _ (3.2) 

By  combining  these  two  equations  we  may  obtain  the  shock  condition 


(3.3) 


In  the  steady  state  situation,  the  reaction  zone  will  be 
seen  to  move  (with  a  constant  velocity  D)  in  a  manner  such  that 
there  is  no  time  variation  of  thermodynamic  variables  in  this 
region.  From  the  equations  of  conservation  for  mass  and  momentum, 

the  detonation  velocity  may  be  expressed  as  D  =  U 


(3.4) 
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In  the  case  of  steady  state  this  formula,  however,  is  valid  not 
only  for  the  pair  (p^,U^)  right  after  the  shock,  but  for  any 
pair  ( p , U)  in  the  reaction  zone.  From  this  one  obtains  a  p,U 
relation  for  positions  in  the  reaction  zone  (the  Rayleight- 
Michelson  line)  [3]  , 

P  =  P0  +  U^T  (U0"U)  .  _ (3.5) 

As  has  been  pointed  out  by  numerous  authors;  see  to  name 
a  few  [  1 ]  —  [ 4 ]  the  three  conservation  equations  (mass,  momentum, 
energy) ,  plus  the  equation  of  state  are  not  sufficient  to 
uniquely  determine  the  detonation  velocity.  Howver,  by  using 
the  Chapman-Jouguet  hypothesis  the  detonation  velocity  is  given 
by 

D  =  u  +  aatS  =  0.  _ (3.6) 

At  a  position  in  the  reaction  zone,  where  the  fraction  of 
nonreacting  molecules  is  3  the  enthalpy  may  be  expressed  as , 

i  =  pU  +  SQ;  _ (3.7) 

if  we  combine  this  with  the  expression  for  the  initial  enthalpy 
(before  detonation)  one  has 

Pu  +  (8*1)  Q  =  -  p (Uq+U)  .  _ (3.8) 
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Combining  these  with  the  p,U  relation  on  the  Rayleigh  line  (and 
assuming  Pq  is  negligibly  small)  we  obtain  p,U  in  terms  of  the 
progress  3,  [4] 


p(S) 


UqTy+TT 


-  Ill 


■-i)Q(i-e) 

T2 - 


(3.9) 


»(«) 


;  I 

Y 


2(y  -l)Q(l-B) 

— 7 - 


It  is  convenient  to  introduce  local  coordinates  in  the  moving 
reaction  zone,  defined  by  a  point  transformation  from  the  (x,t)- 
plane  to  the  (£,t) -plane: 


C  =  tD  -  x, 


X 


(3.10) 


The  variable  £  is  the  distance  from  a  position  x  to  the  shock 
front,  whereas  t  is  the  time  that  has  elapsed  after  the  front 
has  passed  the  point  x.  We  may  express  £  and  t  in  terms  of  8  by 
means  of  a  single  quadrature  as  follows.  From  equation  we  have 


S  =  g(p ,p, 8) 


[ufsT'P(^) '  B  ^  ^  ' 


(3.11) 


consequently  upon  integration  one  has 


t 


r  t 


dx 


dS 

mr 


‘^(0)  , 


(3.12) 


e  -  d^'  -  j  --Itlr3' de  (3) 

J  0  J  1 


(3.13) 


where 

u(B)  =  Ip(0)-PoJ  IUQ-  (6)]  =  p(S)  [  0-  (8)  ].  _ (3.14) 

Elimating  8  between  t  =^(S)  and  £  =  (8),  yields  an  expression 

for  E,  in  terms  of  t,  5  =  ^j  [^>-1  ( t)  ]  . 

In  the  steady-state  it  is  clear  that  the  characteristic 
lines,  I+,  correspond  to  the  straight  lines  8=  constant,  or 
I  :  x  =  tD  +  (6-1) L,  where  L  is  the  reaction  zone  length.  In 
order  to  obtain  I_,  we  return  to  equation  (1.8)  I_:  u(S) 

-  a  (8)  ,  and  compute  a (8)  from 


i2(S)  =  f|§]  =  yp  (8)  v(8)  , 


(3.15) 


1)  (1-8)  Q 


1)  (1-8)0 
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It  is  clear  from  the  above  discussion  that  all  the  "cross- 
characteristics,"  I_,  cut  the  lines  6=constant  at  the  same  angle, 
and  hence  the  I_  form  a  "parallel"  family  of  curves.  From  the 
information  connecting  £ ,  t ,  and  S  we  may  start  at  the  shock 
front,  compute  the  slope  of  I_,  and  then  proceed  by  finite 
differences  across  the  reaction  zone.  The  computation  of  one 
cross-characteristic  will  then  yield  the  entire  family. 

4 .  A  Progressing  Wave  Interpretation: 

In  this  section  we  attempt  to  solve  the  system  of  partial 
differential  equations  (1)  (2)  (3)  (4)  in  terms  of  a  similarity 

variable  [1]  [4]  £  =  xt-a.  We  introduce  this  variable  in  order 

to  reduce  a  system  of  partial  differential  equations  to  a  system 
of  ordinary  differential  equations.*  We  make  the  assumption 
that  the  solutions  u,  p,  p,  have  the  form 


u  *  axt_1U(£) ,, 

p»xkn(£),  _ (4.1) 


and 


P 


a2xk+V2 


P(£)  . 


In  addition  we  introduce. 


and 


S  »  B(£) , 


a2x2t”2  P(£) 

- FT”  JTT TT 


a2x2t‘2 

- FT“ 


E  (£) 


(4.2) 


*We  note,  however,  that  these  solutions  cannot  approach  the  steady- 
state  case  discussed  earlier,  since  here  the  reaction  zone 
continues  to  grow  in  length. 


■r 
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Equations  1.1,  1.2  may  then  be  written  respectively  as 

(U-l)fSJ'  +  KUH  +  mu+SU']  -  0,  _ (4.3) 

U(aU-l)  +  o5U*  (U-l)  +  a  (k+2)  £  +  2^21  =  0,  (4.4) 

ii  ' 

where  primes  indicate  differentiation  with  respect  to  Equation 
1.3  becomes 

QCB '  =  ax2t‘2  2;(,-Qt7^)  E  -  (U-l)Cn'  +  KUO)  ,  (4.5 

Y  £T 

and  then  in  order  for  B  to  be  a  function  of  £  alone  we  must 
have  a  ■  1,  or  ?  *  r.  Equation  1.4  reduces  to 

a£  (U-l) B  '  *  tg ( p  , p , S )  ;  _ (4.6) 


consequently,  for  there  to  exist  a  similarity  solution  (progressing 
wave)  for  our  case  tg(p,p,B)  must  be  expressible  solely  as  a 
function  of  £.  This  amounts  to  a  mathematical  restriction  on 
the  type  of  rate  function  we  may  use;  however,  as  we  shall  see, 
it  does  not  impose  an  important  physical  restriction.  For 
instance,  it  is  customary  to  use  rate  functions  of  the  form  [3] 

[4] 


8  «  cS^p^exp 


(4.7) 


where  eQ  is  an  activation  energy.  In  this  case  equation  (3.6) 


a?  (U-l)  B '  -  tc8m  xk4fll  exp 


Q 

s 


(4.8) 


becomes 
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which  will  be  an  expression  in  £  alone,  if  x  t  =  £  or 
k  =  “-j-  >  1^0.  The  system  (4.1)  -  (4.4)  then  takes  on  the 
form 


(CJ-1)  £Q '  -  ^  fl  +  n[U+£U']  =  0, 
ft(U-l)  [U+£U '  ]  +  £P '  +  (2-  i)P  =  0, 
Q£fl2B'  +  £2P[(U-l)£n'  -  J  Ufl]  =  0, 

and 

2  mi  (i-y )  Eq 

£Z(U-1)B'  *  cBmn£  exp  — - - 2.  . 

£2E(£) 


(4.9) 


This  system  may  be  solved  for  the  first  derivatives  U',  ft', 
P',  B',  as  follows 


U' 


£ (u-i)  {" 


u+cQfl£+1Bm  exp 


un  con*  2Bm 

RTti-IT  “  ~3 - T~  exp 

U  i)  54(u-1)zp 


P'  *  (j  -2)  |  -  cQfi£,+2Bm  exp 


(4.10) 


(4.11) 


(4.12) 


B 


,  _  cb  n 


m  nl  r( 1“Y) i 


£  (0-1) 


exp 


(4.13) 


The  equations  (4.10),  (4.11),  (4.12),  (4.13)  are  a  system  of  four 
ordinary  differential  equations  in  four  dependent  variables. 

The  four  equations  are  essentially  different  as  is  illustrated 
by  solving  for  U',  Q',  P',  B',  and  may  be  solved  by  numerical 


methods . 
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